cutoff <- quantile(ingresos, .95)[[1]]
M = sum(ingresos_censurado>cutoff & ingresos_censurado<top_code)
Te = sum(ingresos_censurado==3)
alpha_est = M/(Te*log(top_code) + sum(log(ingresos_censurado[ingresos_censurado>=cutoff & ingresos_censurado<top_code])) - (M+Te)*log(cutoff))

invertido = alpha_est/(alpha_est-1)
estimated_mean <- top_code*invertido

## Media estimada
estimated_mean

## Media real
real_mean

## ENIGH
## Homogeneizamos los nombres de las columnas de interes
concentradohogar2012 <- read_csv(here("data", "concentradohogar2012.csv"))
colnames(concentradohogar2012)[colnames(concentradohogar2012) == "factor_hog"] <- "factor"

concentradohogar2014 <- read_csv(here("data", "concentradohogar2014.csv"))
colnames(concentradohogar2014)[colnames(concentradohogar2014) == "factor_hog"] <- "factor"

concentradohogar2016 <- read_csv(here("data", "concentradohogar2016.csv"))
concentradohogar2018 <- read_csv(here("data", "concentradohogar2018.csv"))

## Datos anonimos del SAT
sat2012 <- read_delim(here("data", "RUIDO_FINAL_PF_2012.tab"), "\t", escape_double = FALSE, trim_ws = TRUE)

df <- tibble(ing_cor = c(concentradohogar2012$ing_cor,
                         concentradohogar2014$ing_cor,
                         concentradohogar2016$ing_cor,
                         concentradohogar2018$ing_cor),
             factor = c(concentradohogar2012$factor,
                        concentradohogar2014$factor,
                        concentradohogar2016$factor,
                        concentradohogar2018$factor),
             año = c(rep("2012", length(concentradohogar2012$ing_cor)),
                     rep("2014", length(concentradohogar2014$ing_cor)),
                     rep("2016", length(concentradohogar2016$ing_cor)),
                     rep("2018", length(concentradohogar2018$ing_cor)))
              )

histograma <- ggplot(df, aes(x = ing_cor, weight = factor)) + 
  geom_histogram(bins = 60) + 
  facet_wrap(~año) +
  scale_x_continuous(labels = scales::comma) +
  scale_y_continuous(labels = scales::comma) +
  xlab("Ingreso Corriente") +  ylab("Frecuencia") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

histograma_zoomed <- ggplot(df, aes(x = ing_cor, weight = factor)) + 
  geom_histogram(bins = 600) + 
  facet_wrap(~año) +
  scale_x_continuous(labels = scales::comma) +
  scale_y_continuous(labels = scales::comma) +
  coord_cartesian(xlim=c(0,3000000)) +
  xlab("Ingreso Corriente") +  ylab("Frecuencia") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

histograma_log <- ggplot(df, aes(x = ing_cor, weight = factor)) + 
  geom_histogram(bins = 30) + 
  facet_wrap(~año) +
  scale_x_continuous(trans = "log", labels = scales::comma) +
  scale_y_continuous(labels = scales::comma) +
  xlab("Ingreso Corriente (log)") +  ylab("Frecuencia") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

ingreso_medio_percentil <- function(datos, n_percentil) {
  # datos son las base de datos del censo, n_percentil es el numero de divisiones (e.g.: 10 retorna deciles, 100 retorna percentiles)
  
  ## Datos sorteados por ingresos trimestales corrientes más indice de los hogares representados observados
  sorted_ing <- datos %>% select(ing_cor, factor) %>%
    arrange(ing_cor) %>%
    mutate(cumulative_obs = cumsum(factor))
  ## Total de hogares representados (ver definicion de factor en documentacion de la ENIGH)
  total_hogares <- sum(datos$factor)
  
  ## Tamaño por decil truncado para no incluir decimales
  tam_decile <- trunc(total_hogares/n_percentil)
  
  ## Primeros n-1 deciles
  lim_inf <- 1
  lim_sup <- tam_decile
  ing_med_de <- c()
  for (i in 1:(n_percentil-1)) {
    ing_med_de[i] <- sorted_ing %>% filter(cumulative_obs >= lim_inf, cumulative_obs <= lim_sup) %>%
      summarise(mean = sum(ing_cor*factor)/tam_decile) %>% 
      pull()
    
    ## Limites del decil
    lim_inf <- lim_inf + tam_decile
    lim_sup <- lim_sup + tam_decile
  }
  
  ## Decimo decil aparte para incluir el resto de observaciones
  ing_med_de[n_percentil] <- sorted_ing %>% filter(cumulative_obs >= lim_inf) %>%
    summarise(mean = sum(ing_cor*factor)/tam_decile) %>% 
    pull()
  
  return(ing_med_de)
}

n_percentil <- 10
df <- tibble("2012" = ingreso_medio_percentil(concentradohogar2012, n_percentil),
             "2014" = ingreso_medio_percentil(concentradohogar2014, n_percentil),
             "2016" = ingreso_medio_percentil(concentradohogar2016, n_percentil),
             "2018" = ingreso_medio_percentil(concentradohogar2018, n_percentil),
             index = 1:n_percentil)

df <- df %>% gather(key = "Año", value = "Ingreso_medio", -index)

ingreso_medio_decil <- ggplot(df , aes(x = index, y = Ingreso_medio))+
  geom_bar(stat="identity")+
  facet_wrap(~ Año) +
  scale_y_continuous(labels = scales::comma, breaks = seq(0,175000,25000))+
  scale_x_continuous(breaks = 1:10, labels = scales::ordinal)+
  xlab("Decil") + ylab("Ingreso medio") +
  labs(title="Ingreso medio por deciles") +
  theme_bw()

ingreso_medio_decil

ingreso_acumulado_percentil <- function(datos, n_percentil){
  # Retorna el porcentaje del ingreso total que acumula cada percentil
  # datos son las base de datos del censo, n_percentil es el numero de divisiones (e.g.: 10 retorna deciles, 100 retorna percentiles)
  
  ## Datos sorteados por ingresos trimestales corrientes más indice de los hogares representados observados
  sorted_ing <- datos %>% select(ing_cor, factor) %>%
    arrange(ing_cor) %>%
    mutate(cumulative_obs = cumsum(factor))
  
  ## Total de hogares representados (ver definicion de factor en documentacion de la ENIGH)
  total_hogares <- sum(datos$factor)
  
  ## Muestras por por percentil truncado para no incluir decimales
  tam_decile <- trunc(total_hogares/n_percentil)
  
  
  ## Primeros n-1 deciles
  lim_inf <- 1
  lim_sup <- tam_decile
  ing_total <- c()
  for (i in 1:(n_percentil-1)) {
    ing_total[i] <- sorted_ing %>% filter(cumulative_obs >= lim_inf, cumulative_obs <= lim_sup) %>%
      summarise(mean = sum(ing_cor*factor)) %>% 
      pull()
    
    ## Limites del decil
    lim_inf <- lim_inf + tam_decile
    lim_sup <- lim_sup + tam_decile
  }
  
  ## Decimo decil aparte para incluir el resto de observaciones
  ing_total[n_percentil] <- sorted_ing %>% filter(cumulative_obs >= lim_inf) %>%
    summarise(mean = sum(ing_cor*factor)) %>% 
    pull()
  
  ing_total <- ing_total/sum(datos$ing_cor*datos$factor)
  
  return(ing_total)
}

n_percentil <- 10
df <- tibble("2012" = ingreso_acumulado_percentil(concentradohogar2012, n_percentil),
             "2014" = ingreso_acumulado_percentil(concentradohogar2014, n_percentil),
             "2016" = ingreso_acumulado_percentil(concentradohogar2016, n_percentil),
             "2018" = ingreso_acumulado_percentil(concentradohogar2018, n_percentil),
             index = fct_rev(as.factor(1:n_percentil))) 

df <- df %>% gather(key = "Año", value = "Ingreso", -index)

ingreso_acu_decil <- ggplot(df, aes(fill = index, x = Año, y = Ingreso)) +
  geom_bar(position="stack", stat="identity") +
  ylab("Porcentaje de Ingreso acumulado") + xlab("Año") +
  scale_fill_viridis(discrete = T, name = "Decil") +
  scale_y_continuous(labels = scales::percent) +
  theme_bw() +
  labs(title= "Ingreso acumulado por Decil")+
  theme(legend.background = element_rect(fill = "transparent", colour = "transparent"))

ingreso_acu_decil

n_percentil <- 1000
df <- tibble("2012" = ingreso_acumulado_percentil(concentradohogar2012, n_percentil),
             "2014" = ingreso_acumulado_percentil(concentradohogar2014, n_percentil),
             "2016" = ingreso_acumulado_percentil(concentradohogar2016, n_percentil),
             "2018" = ingreso_acumulado_percentil(concentradohogar2018, n_percentil),
             index = 1:n_percentil)

df <- df %>% gather(key = "Año", value = "Ingreso", -index)

ingreso_acu_per <- ggplot(df, aes(x = index, y = Ingreso)) +
  geom_line(aes(color = Año), size = 1) +
  scale_x_continuous(labels = scales::ordinal) +
  scale_y_continuous(labels = scales::percent) +
  coord_cartesian(xlim=c(990,1000)) +
  ylab("Porcentaje de Ingreso acumulado") + xlab("Percentil (0.1%)") +
  labs(title="Ingreso acumulado por cada 0.1% de la poblacion") +
  theme_bw() +
  theme(legend.position = c(.1,.8), legend.background = element_rect(fill = "transparent", colour = "transparent"))

ingreso_acu_per

gastan_de_mas <- function(datos){
  total_hogares <- sum(datos$factor)
  datos %>% 
    transmute(mas_gasto = (gasto_mon > ing_cor)*factor) %>%
    summarise(proporcion = sum(mas_gasto)/total_hogares) %>%
    pull()
}

# Porcentaje de personas que gastan de más
over_spenders <- c(gastan_de_mas(concentradohogar2012),
                   gastan_de_mas(concentradohogar2014),
                   gastan_de_mas(concentradohogar2016),
                   gastan_de_mas(concentradohogar2018))

over_sp <- ggplot()+
  geom_bar(aes(x=c("2012", "2014", "2016", "2018"), y=over_spenders), stat = "identity") +
  scale_y_continuous(labels = scales::percent) +
  xlab("Año") + ylab("")+
  labs(title="Porcentaje personas que gastan más de lo que ingresan")+
  theme_bw()

over_sp 

breaks <- c(1,2,3,4,5,6,8,10, 15,20,30,50,100,200,500,1000)*1000000


means <- c()
for (i in 1:(length(breaks)-1)) {
  means[i]<- mean(sat2012$I_DEC_TIAONCT1_AA[sat2012$I_DEC_TIAONCT1_AA >= breaks[i]
                                            & 
                                            sat2012$I_DEC_TIAONCT1_AA<breaks[i+1]], 
                  na.rm = T)
}

top_sat <- ggplot()+
  geom_line(aes(x=1:(length(breaks)-1), y=means), color = "light blue", size = 1) +
  scale_x_continuous(breaks = 1:(length(breaks)-1), labels = c("1-2 M", "2-3 M", "3-4 M", "4-5 M", "5-6 M", "6-8 M", "8-10 M", "10-15 M", "15-20 M", "20-30 M", "30-50 M", "50-100 M", "100-200 M", "200-500 M", "500-1000 M")) +
  scale_y_continuous(labels = scales::comma, breaks=seq(0,600000000,100000000)) +
  labs(title = "Ingreso medio segun ingresos reportados al SAT") +
  xlab("Ingreso medio") + ylab("Rango de ingreso")+
  theme_bw() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

top_sat

mean_log <- function(x) {
  exp(x[1]+(x[2]^2)/2)
}

pll_log <- function(x){
  -sum(dlnorm(concentradohogar2012$ing_cor+.01, x[1], x[2], log = TRUE)*concentradohogar2012$factor)
}

par_log_NR<- optim(par = c(10, 1), pll_log, control=list(maxit=1000))$par


par_log_R <- solnp(pars = c(11, 2), # Valores inciales
               pll_log, # función a optimizar
               eqfun=mean_log, # función restricción
               eqB=92733.62,   # Valor a igualar
               LB=c(0,0), # Soporte mínimo para los parametros
               UB=c(20,5))$pars # Soporte máximo para los parametros

## Media muestral
weighted.mean(concentradohogar2012$ing_cor, concentradohogar2012$factor)

## Media No restringida
mean_log(par_log_NR)

## Media restringida
mean_log(par_log_R)

mean_gb2 <- function(x) {
  moment.gb2(1, x[1], x[2], x[3], x[4])
}

pll_gb2 <- function(x){
  -loglp.gb2(concentradohogar2012$ing_cor+.01, x[1], x[2],x[3], x[4], w=concentradohogar2012$factor)
}

par_gb2_NR<- optim(par = c(1.307071, 20594.790235, 2.341514, 1.843004), pll_gb2, method = "BFGS",control=list(maxit=1000))$par

par_gb2_R <- solnp(pars = c(1.5, 30000, 1, 1), # Valores inciales
               pll_gb2, # función a optimizar
               eqfun=mean_gb2, # función restricción
               eqB=92733.62,   # Valor a igualar
               LB=c(0,0,0,0), # Soporte mínimo para los parametros
               UB=c(5,100000,5,5))$pars # Soporte máximo para los parametros

## Media muestral
weighted.mean(concentradohogar2012$ing_cor, concentradohogar2012$factor)

## Media No restringida
mean_gb2(par_gb2_NR)

## Media restringida
mean_gb2(par_gb2_R)

distribucion_empirica <- function(x, datos, factor){
  # x es el valor donde se quiere evaluar la función de distribución empírica
  # datos es el conjunto de datos en formato de vector que se quiere evaluar
  df <- data.frame(datos = datos, factor = factor, cumu = cumsum(factor))
  F_x <- c()
  for (i in 1:length(x)) {
    F_x[i] <- sum(df[df$datos<x[i], "factor"])/df[length(datos), "cumu"]
  }
  return(F_x)
}
distri_empi_restri <- ggplot() +
  geom_function(aes(colour = "Empirica"), 
                fun = distribucion_empirica, 
                args = list(datos = concentradohogar2012$ing_cor, 
                            factor = concentradohogar2012$factor), 
                size = 1) +
  geom_function(aes(colour = "GB2"), 
                fun = pgb2, 
                args = list(shape1 = par_gb2_R[1],
                            scale = par_gb2_R[2],
                            shape2 = par_gb2_R[3],
                            shape3 = par_gb2_R[4]),
                size = 1) +
  geom_function(aes(colour = "Log-Normal"), 
                fun = plnorm, 
                args = list(meanlog = par_log_R[1],
                            sdlog = par_log_R[2]), 
                size = 1) +
  scale_x_continuous(trans = 'log', limits = c(600, 4000000), breaks = 10^(seq(2,6,1)), labels = scales::comma) +
  scale_y_continuous(breaks = seq(0,1,.1)) +
  scale_colour_manual(name="Distribución",values=c("Empirica"="#3591d1", "GB2"="#f04546", "Log-Normal"="#f045d6")) +
  ylab("Probabilidad") + xlab("Ingreso (Escala log)") +
  labs(title = "Distribución empírica contra estimadas con restricción") +
  theme_bw() +
  theme(legend.position = c(.1,.8), legend.background = element_rect(fill = "transparent", colour = "transparent"))

distri_empi_restri

p1 <- ggplot() +
  geom_function(aes(colour = "Log-Normal Restringida"), 
                fun = plnorm, 
                args = list(meanlog = par_log_R[1],
                            sdlog = par_log_R[2]), 
                size = 1) +
  geom_function(aes(colour = "Log-Normal No Restringida"), 
                fun = plnorm, 
                args = list(meanlog = par_log_NR[1],
                            sdlog = par_log_NR[2]), 
                size = 1) +
  scale_x_continuous(trans = 'log', limits = c(600, 4000000), breaks = 10^(seq(2,6,1)), labels = scales::comma) +
  scale_y_continuous(breaks = seq(0,1,.1)) +
  scale_colour_manual(name="Distribución",values=c("Log-Normal Restringida"="#f04546", "Log-Normal No Restringida"="#3591d1")) +
  ylab("Probabilidad") + xlab("Ingreso (Escala log)") +
  theme_bw() +
  theme(legend.position = c(.25,.8), legend.background = element_rect(fill = "transparent", colour = "transparent"))

p2 <- ggplot() +
  geom_function(aes(colour = "GB2 Restringida"), 
                fun = pgb2, 
                args = list(shape1 = par_gb2_R[1],
                            scale = par_gb2_R[2],
                            shape2 = par_gb2_R[3],
                            shape3 = par_gb2_R[4]),
                size = 1) +
  geom_function(aes(colour = "GB2 No Restringida"), 
                fun = pgb2, 
                args = list(shape1 = par_gb2_NR[1],
                            scale = par_gb2_NR[2],
                            shape2 = par_gb2_NR[3],
                            shape3 = par_gb2_NR[4]),
                size = 1) +
  scale_x_continuous(trans = 'log', limits = c(600, 4000000), breaks = 10^(seq(2,6,1)), labels = scales::comma) +
  scale_y_continuous(breaks = seq(0,1,.1)) +
  scale_colour_manual(name="Distribución",values=c("GB2 Restringida"="#f04546", "GB2 No Restringida"="#3591d1")) +
  ylab("Probabilidad") + xlab("Ingreso (Escala log)") +
  theme_bw() +
  theme(legend.position = c(.2,.8), legend.background = element_rect(fill = "transparent", colour = "transparent"))

grid.arrange(p1,p2,ncol=2)

## Calculamos el acumulado por decil mediante simulación
x <- rlnorm(100000, par_log_R[1], par_log_R[2])
deciles <- quantile(x, probs = seq(0,1,.1))
total<-sum(x)

medd <- c()
for (i in 1:10) {
  medd[i] <- sum(x[x > deciles[i] & x <= deciles[i+1]])/total
}

## Construcción de la tabla para alimentar el ggplot
n_percentil <- 10
df <- tibble("ENIGH 2012" = ingreso_acumulado_percentil(concentradohogar2012, n_percentil),
             "Log-Normal" = medd,
             index = fct_rev(as.factor(1:n_percentil))) 

df <- df %>% gather(key = "Fuente", value = "Ingreso", -index)

## Construcciónd de la gráfica
ingreso_acu_decil_ENIGH_Log <- ggplot(df, aes(fill = index, x = Fuente, y = Ingreso)) +
  geom_bar(position="stack", stat="identity", width=0.7) +
  ylab("Porcentaje de Ingreso acumulado") + xlab("Modelo") +
  scale_fill_viridis(discrete = T, name = "Decil") +
  scale_y_continuous(labels = scales::percent) +
  theme_bw() +
  theme(legend.background = element_rect(fill = "transparent", colour = "transparent"))

ingreso_acu_decil_ENIGH_Log